clear;
clc;

cd('specify directory here');
load('EdEx_matlab');

global Y T1 T2 SM SR Draws M Xc or orm Gmax Gmin ll 
options=optimoptions('fminunc','Algorithm','quasi-newton','Display','iter','PlotFcn',{@optimplotx,@optimplotfval,@optimplotfirstorderopt});

rng(10010);
Drawso=normrnd(0,1,6063,500);
M=500;

%% WHITES
I=(cbyrace6==1)+0;

Y=Yo(I==1);
T1=T1o(I==1);
T2=T2o(I==1);
SM=SMo(I==1);
SR=SRo(I==1);
Draws = Drawso(I==1,:);
Xc=Xco(I==1);
or = repmat((cbyerace6(I==1)==0)+0,1,M);
orm = repmat((cbymrace6(I==1)==0)+0,1,M);

theta2=[0;0;1;1;1;1;0;1;1;0.6;0.6;1;1;0.2;0.2;1;0.5;0.3;0.5;0.1;0.1;0.2;0.2;0.1;0.1;1;1];
eps=1e-5;
no = 0.1;
while no>eps
    theta = theta2;
    [theta2, FVAL]=fminunc(@(param)LLcorr_probit(param),theta,options);
    no=norm(theta2-theta);
end

PARAMwhite = theta2;
varPARAMwhite = ERRORS_corr(PARAM);
sewhite = sqrt(diag(varPARAM));

%% BLACKS
I=(cbyrace3==1)+0;

Y=Yo(I==1);
T1=T1o(I==1);
T2=T2o(I==1);
SM=SMo(I==1);
SR=SRo(I==1);
Draws = Drawso(I==1,:);
Xc=Xco(I==1);
or = repmat((cbyerace3(I==1)==0)+0,1,M);
orm = repmat((cbymrace3(I==1)==0)+0,1,M);

theta2=[0;0;1;1;1;1;0;1;1;0.6;0.6;1;1;0.2;0.2;1;0.5;0.3;0.5;0.1;0.1;0.2;0.2;0.1;0.1;1;1];
eps=1e-5;
no = 0.1;
while no>eps
    theta = theta2;
    [theta2, FVAL]=fminunc(@(param)LLcorr_probit(param),theta,options);
    no=norm(theta2-theta);
end

PARAMblack = theta2;
varPARAMblack = ERRORS_corr(PARAMblack);
seblack = sqrt(diag(varPARAMblack));

